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Abstract 

A 3-D numerical model of comet nuclei is presented. An implicit numerical scheme 
was developed for the thermal evolution of a spherical nucleus composed of a mix- 
ture of ice and dust. The model was tested against analytical solutions, simplified 
numerical solutions, and 1-D thermal evolution codes. The 3-D code was applied to 
comet 67P/Churyumov-Gerasimenko; surface temperature maps and the internal 
thermal structure was obtained as function of depth, longitude and hour angle. The 
effect of the spin axis tilt on the surface temperature distribution was studied in 
detail. It was found that for small tilt angles, relatively low temperatures may pre- 
vail on near-pole areas, despite lateral heat conduction. A high-resolution run for a 
comet model of 67P/Churyumov-Gerasimenko with low tilt angle, allowing for crys- 
tallization of amorphous ice, showed that the amorphous /crystalline ice boundary 
varies significantly with depth as a function of cometary latitude. 
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1 Introduction 



Our knowledge of the nature of comets, their structure, composition, dynam- 
ical history, and modes of activity has greatly improved over the last twenty 
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years. The one event that constituted a mile-stone in the study of comets was 
the space mission Giotto, which provided the first close-up picture of a comet's 
heart - its nucleus. Several other space m issions to comets foll owed, crowned 



by the very recent Deep Impact mission (lA'Hearn et all 120051 ) . aimed at re- 
vealing properties of the nucleus interior. As observational data on comets has 
been accumulating, theoretical studies and modeling have started to develop, 
in order to interpret, predict and explain new findings and discoveries. 



Comets consist of frozen gases and dust. The dust consists of silicates and or- 
ganic materials, while the frozen gases are mostly water and a few constituents 
more volatile than water. Solar heat, not reflected or radiated from the sur- 
face of the nucleus is primarily used to evaporate ices from exposed areas. The 
remainder of the heat penetrates into the nucleus, where it can cause numer- 
ous physico-chemical reaction processes, which eventually affect measurable 
properties at the nucleus surface. Over the past two decades, numerical stud- 
ies have attempted to integrate some of the more influential processes with 
a numerical representation of a comet nucleu s, in order to better understand 



the rather unpredictable behavior of comets (IPrialnik et all 120041 ) 



The most obvious property of comet nuclei is the lack of spherical symmetry. 
Besides the fact that self-gravity, which dictates the spherical shape of larger 
celestial bodies, is negligible in comets, the main energy source — solar radi- 
ation — is unevenly distributed over the nucleus surface, not only because of 
daily variation and spin axis inclination, but also because of the large orbital 
variations, characteristic to comets, but not to other bodies of the solar sys- 
tem. Nevertheless, theoretical models of nuclei have so far assumed spherical 
symmetry, or have circumvented this assumption by crude approximations. 

Modeling began with very simple pictures of the nucleus. A few physical pro- 



for the comet nucleus behavior ( 


Cowan and Ahearn. 


1979; 


Fanale and Sal vail. 


1984; 


Herman and Podolak. 


1985|; 


Prialnik and Bar- Nun. 1987h. However, these 



models were able to give only average results and explain only the general be- 
havior of comets. The models were limited to a single dimension, mostly based 
on spherical symmetry (the fast rotator approximation). They simulated heat 
flow in a mixture of ice and rock (dust) with uniform boundary conditions 
for the solar flux, radiative emission and substance sublimation. Additional 
physico-chemical processes, e.g., amorphous ice crystallization, surface ero- 
sion etc., were added later on to produce a more elaborate physical picture. 
As these models did not take into consideration spatial variations in surface 
temperature, they provided a fair approximation only for the internal part of 
the nucleus, where these variations fade. 



For more realis tic surface results, s l ow rotator approache s were adopted. A 
"1.5-D" model (IBenkhoff and Boicd . Il996l ; iBenkhofil . Il999l ) was suggested, in 



2 



which, a 1-D model was applied with boundary conditions of an equatorial 
point on a rotating nucleus. This way, diurnal boundary conditions differ- 
ences were taken into considerations, and an upper limit for production rate 
was achieved when the sub-solar (high noon) flux was adopted for the entire 
surface of the sunlit hemisphere. This model, as its predecessor, does not cal- 
culate lateral flow. A sem i-lateral approach was taken for the 2.5-D model 
( lEnzian et all 119971 . Il999l ) where a meridian flow was calculated, and bound- 
ary conditions were taken as in the 1.5-D model. 



A dif f erent approach is th e quasi-3-D procedure ([Gutierrez et all l2000l ; I Julian et al. 



2000l ; ICohen et all 120031 ) . In this approach, every point on the nucleus surface 



is calculated (with the appropriate local boundary condition) with radial flow 
only, so calculations are similar to the 1.5-D method, but boundary condi- 
tions are of the total sphere, just as with a 3-D approach. Still, each point (or 
element) of the surface evolves independently of the others and only radial 
conduction is considered. 



Now that modeling has made significant progress towards understanding the 
general structure and behavior of comets, more sophisticated and realistic 
pictures of the nucleus are required. 



The purpose of the present research project is to develop a fully 3-D model of a 
comet nucleus. As in the quasi-3-D, boundary conditions will be taken for the 
entire sphere, but now, meridional and azimuthal heat fluxes are calculated 
as well. Gas production and flow in the interior of the comet are not included 
in the present model. In order to compare production rates of volatiles with 
the observed ones at large distances from the nu c leus, when the surface is 



not resolved (in fact, not even seen) (IBiver et all Il999l ). a sum of all local 



production rates can be calculated. However, the model will provide a detailed 
description of the surface activity as well. 



In the 3-D model, in contrast to the older models, all the grid elements need 
to be solved simultaneously, so the computational load increases dramatically. 
This may limit the spatial resolution of the comet if we require a solution in 
a reasonable time. However, as modern computers advance rapidly, the grid 
density may increase for the same calculation time. This model may be scaled 
for super-computers to allow even denser grid resolution. 



2 Model description 



Our purpose is to calculate numerically a variety of physical processes that 
are believed to take place inside a comet nucleus. In order to do so, physical 
equations and models should be adjusted for discrete calculations. The basic 
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adjustment is dividing the nucleus into elements via a grid, and simplifying 
the problem by assuming a "homogeneous lumped system" approximation for 
each element. 

2.1 3-D heat conduction equation 

Heat conduction is modeled by Fourier's law 

F = -KVT, (1) 

where F is the heat flux, K is the conduction coefficient and T is the temper- 
ature. 

In 3-D spherical coordinates, equation (1) will take the form 

\or r ad r sin V o<p J 

where H z is the Hertz factor, a correction factor (< 1) for the heat conductivity 
of porous substances that accounts for the reduced contact area between solids 
as compared to the cross-sectional area. In order to calculate the effective 
conduction coefficient K for a "homogeneous lumped system" we considered 
all substances distributed evenly through the system, and the contribution 
of each substance to the final value is taken to be proportional to its mass 
fraction X a . Thus 

K = Y J X a K a , (3) 

a 

where K n is the conduction coefficient of substance a. 



When all fluxes through a "lumped system" are combined we get the heat 
equation 

-V-F + q = pC p ^, (4) 

where q is the rate of energy release by an internal heat source within the 
"lumped system" , p is the bulk density and C p is the heat capacity. 



2.1.1 Amorphous to crystalline ice transition 

The crystallization process of amorphous ice is spontaneous but greatly af- 
fected by temperature. Its rate is given by (ISchmitt et all Il989l ). based on 
laboratory experiments 



A 



1.05 x 10 13 < 



-5370/T 



(5) 
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The amorphous ice mass fraction X A changes with time at a rate 

X A = -X(T)X A (6) 

This process, which changes the ice structure, also changes its physical prop- 
erties, such as density and thermal conductivity. In addition, crystallization 
is an exothermic process that can cause chain-reaction, that is, propagate by 
feeding on its own energy. The rate of heat generated by crystallization is 
given by 

q = X(T)pX A H a (7) 



where 7i a is the energy released per unit mass flGhormleyl . Il968l ). 



2.2 Boundary conditions 



Boundary conditions at the surface consist of balance of heat fluxes into the 
nucleus from external sources and out of it into space. Three such sources are 
considered: 



'1) Solar radiation flux, given by 



(1 -^4^ C0S? ' 



where A is the surface albedo, L Q is the solar luminosity, du is the he- 
liocentric distance and £ is the local solar zenith angle. Here cos£ is 
calculated for surface coordinates 9 and <fi according to 

cos £ = sin 9 S sin 9 cos ((f) a + </>)+ cos 9 S cos 9, (9) 

where 9 S and cf> s are the sub-solar coordinates. 

(2) Thermal radiation loss, given by 

eoT\ (10) 

where e is the emissivity and a is the Stefan-Boltzmann constant. 

(3) Latent heat absorbed per unit time by sublimation of volatiles, given by 



where m is the molecular weight, H is the latent heat of sublimation and 
V v is the vapor pressure, given by the Clausius-Clapeyron equation 

V ae -b/T r 12 \ 
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Fig. 1. Logical representation of the numerical grid, conforming to spherical coor- 
dinates convention. The O coordinate is divided into 'J' equal intervals (AO), the 
<5 coordinate is divided into , K' equal intervals (A<1?) and the radial coordinate is 
divided with geometric progression into intervals (ARi). The equations of the 
numerical scheme are based on this representation. 

2.3 Numerical scheme 

The heat equation (4) is a non-linear second-order partial differential equa- 
tion, which has no analytical solution and has to be solved numerically. The 
differential equation is transformed into a set of difference equations, and ap- 
plied to a finite grid, where the infinitesimals 5R,5Q,5& and St become the 
grid steps AR, A0, A$ and At, respectively. 

The grid has two representations, the logical (geometric) one and the com- 
puterized (algebraic) one. The logical grid, shown in Fig. 1, is the spatial 
representation of the difference equations, and is a 3-D grid in spherical coor- 
dinates. The grid itself is a 3-D mesh with I, J and K divisions for R (radial 
distance), (co-latitude) and <3> (azimuth). In order to better study surface 
phenomena, AR is taken to be decremental geometrically, and is uniquely de- 
termined by the radius of the nucleus, the number of divisions ('/') and the 
surface layer's thickness. Dimensions G and $ are divided into equal intervals. 

The computerized grid is a vector with n = I-J-K elements. The elements are 
organized along the vector, as shown in Fig. 2, so that radius values will vary 
the slowest, and the co-latitude will vary the fastest. The relation between the 
logical and the computerized grid is expressed by: 

v = j + (k-l) J+(i-l) JK (13) 

v is the position in the computer grid vector, while i,j and k are the indexes 
for R,Q and $ and J and K are number of divisions for O and $. This 
configuration forms a diagonal equation matrix, which is easier to solve. 
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Fig. 2. Computerized (vector) representation of the numerical grid. This represen- 
tation is the order of the grid's elements in the solution vector of size n = I J K. 
The relation between the logical and vector grids is given by equation (13). In this 
representation R varies the slowest through the vector and varies the fastest. 
The equations matrix, according to this representation, will be diagonal (albeit not 
tight) and easier so solve. 



2-4 Difference equations 



For greater stability, and in order to remove time step constraints, a fully im- 
plicit iterative scheme is chosen, although the formulation is more complicated 
and has a higher computational load. The difference equation for each volume 
element has the form 

Y, F^ n) AS x + q^AV = — AV (14) 

where x is the direction (i + ,i~ ,j + ,j~ ,k + and k~), n is the iteration number, F x 
is the heat flux vector component in the x direction, AS X is the element's re- 
side area, q n is produced heat inside the element, AV is the element's volume, 
At is the time step, and U is the thermal energy. We note that the difference 
scheme is conservative, that is, upon integration over the entire volume, all 
fluxes cancel out except the flux crossing the nucleus surface. In other words, 
Gauss's theorem is satisfied by the difference equations on the discrete grid. 
The sum on the LHS of equation (14) runs over all sides of a volume element: 
6 for most elements, 5 for those ending on the = or n axis, and 4, for 
those ending at the center R = 0. 

The equations are linearized and solved for AT rather than for T itself; as AT 
reduces over iterations, it serves as convergence criterion. 
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3 Tests of the 3-D model 



3. 1 Analytical and numerical tests 



Several comparative tests were carried out to ensure that heat flow is correctly 
solved. First, two analytical solutions for the heat transport equation on a 
sphere were chosen. Since all terms in our surface boundary condition (solar 
radiation, sublimation and thermal emission) are non-linear, and hence do not 
allow analytical solutions to the heat equation, we chose a linear convective 
boundary condition. The numerical model was modified to the same boundary 
condition, and kept spherically symmetric in order to become comparable to 
the analytical soluti on. The first analytical mo del (hereafter, C&J) was taken 
from the literature (jCarslaw and Jaegerl . Il959l ): 



T(r,t) 



2hT n 



E 



sin (A n a) (A n a) 2 + (1 — ah)' 



X„ defined as: 



XI [XI - ah (1- ah)] 
tan (X n a) 



sin (X n r) e 



-XiDt 



(15) 



(16) 



1 — ah 

where T is the initial temperature, a is the sphere's radius, h is the convection 
constant, D is the conduction constant and t is the time step; X n is the series 
solution of Eq.(16). This solution is based on simplifying assumptions. 



A second, accurate analytical solution was derived: 



T(r,t) 



2T ^ sin (X n a) — (X n a) cos (X n a) 



E 



X^a — X n sin (X n a) cos (X n a 



sin (X n r) e 



-XiDt 



(17) 



which is independent of h, and therefore not limited by boundary condition 
values. 

A different test model was obtained by solving the heat equation numerically 
by a 1-D spherical explicit scheme: 



T 



(n+l) rp{n) 



TP + 2D 



At 



rp{n) _ rp(n) rp{n) _ rp(n) 

- L i—1 -*t qi _|_ i+1 i c< 

Ar { + Arj_x 1 Ar { + Ar i+1 t+1 



with boundary condition (i = I): 



(n+l) 



T / (n) + 2D 



At 



rp{n) rp{n) 1 



(19) 
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Fig. 3. Different solutions for a cooling sphere with uniform initial conditions, and 
linear convection boundary condition, at 50, 1000 & 10000 seconds. All solutions: 
Carslaw & Jaegar, analytical development, numerical 1-D and the 3-D (Komet3D) 
models converge to the same result (left). The absolute difference between all models 
to the 3-D model (right) reveals perfect agreement between the numerical models, 
and up to 5% error for the analytical solutions. 

where 



-7T 



(Ri + dRi) 3 — R\ 



and 



Si = 4ttW 



(20) 



(21) 



The tests were performed on a sphere with R = 1000 cm. All other parameters 
were chosen artificially (to allow the C&J model to converge) and do not rep- 
resent a comet nucleus. The sphere was initially set at a uniform temperature 
of 60°K and was then allowed to cool. Temperature profiles as function of 
depth are shown in Fig. 3 at three different times for all 4 solutions. Clearly, 
the 3-D model agrees very well with both analytical models, as well as with the 
1-D numerical model at all times. A more detailed comparison of the model 
with the other models is obtained by plotting the differences between them, 
as shown in the right panel of Fig. 3. The maximal difference between the 
3-D model and the analytical models is less than 5%, which is an excellent 
agreement, keeping in mind that the lowest radial grid resolution is 9%; the 
difference between the two numerical models is negligibly small. 



3.2 Comparison with other models for comet J^6P '/Wirtanen 



Comet 46P /Wirtanen was intensively studied while it was the target of space 
mission Rosetta. In particular, several different and independent numerical 
codes were applied to this comet and the results were carefully analyzed and 
compared (Huebner et al. 1999). These models, allowing for differences among 
them, can be considered a reliable test for the present model. We therefore 
used our code adopting comet 46P/Wirtanen's parameters (a=3. 093738 AU, 
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Fig. 4. Orbital temperature variation of an equatorial model element with comet 
Wirtanen's orbital parameters, for 5 revolutions around the sun. Colors distinguish 
between revolutions. Line thickness indicates temperature differences between day 
and night sides. 



e=0. 6578222, R=600m, and a spin period of 24 hr, as used in other model 
calculations), and calculated the thermal evolution for 5 revolution around 
the sun, assuming the spin axis to be perpendicular to the orbital plane (zero 
tilt angle). The nucleus composition consisted of 50% crystalline H 2 ice and 
50% dust, by mass. No Hertz factor was used (H z = 1). In Fig. 4, the temper- 
ature of a point on the equator is shown versus time and heliocentric distance. 
The line thickness represents daily temperature variations. As expected, the 
perihelion (daytime) temperature is the same for all orbits, since heat ex- 
change with the interior is negligible in this case, and the surface temperature 
is controlled by sublimation. The aphelion temperature increases slightly with 
each period; again, this is to be expected, as the nucleus surface layer ac- 
cumulates heat and needs several revolutions to stabilize. The results are in 
very good agreement with Huebner et al. (loc. cit.), for a model with sim- 
ilar parameters. We also compared our results with th e quasi-3-D model of 
comet 46P/Wirtanen computed by ( Cohen et al. . 20031 ) and found very good 
matching of the temperature map. 



Fig. 5 illustrates boundary fluxes for that same point on the equator. Positive 
fluxes are directed inward, and negative fluxes outward. At perihelion, the 
solar flux is mostly absorbed by sublimation, while at aphelion the various 
terms are comparable. The residual flux, which represent heat conducted into 
or out of the nucleus, shows that on the day side, there is a net inward flux, 
and on the night side, a net outward flux. 
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Time [days] 



Fig. 5. Surface flux balance for equatorial element, with comet Wirtanen's orbital 
parameters (same as in Fig. 4). Positive values are inward fluxes, and negative are 
outwards fluxes (relative to the nucleus center). The red line represents total flux. 

4 Applications of the 3-D code 



4-1 Models of Comet 67P/ Churyumov-Gerasimenko 



The actual target of the Rosetta mission is now comet 67P/ Churyumov- 
Gerasimenko (hereafter, 67P/C-G). We have used our code to model the 
evolution of a spherical nucleus that has the characteristics of this comet, 



a=3 5029497 AU, e=0.6319359, R=1980m, and a spin period o f 12.6 hr (ILamv et al. 



20041 ). in particular a spin axis tilt of 40°-45° f lChesleyl . 120041 ) . We considered 
a nucleus containing 50% water crystalline ice and 50% dust by mass, and 
calculated the thermal evolution for 5 orbital revolutions. Figures 6 and 7 
illustrate the daily and orbital skin depths, by comparing several tempera- 
ture maps for different depths and orbital positions. Analytical diurnal and 
orbital skin-depths were calculated to be 11 cm and 8 m respectively, in good 
agreement with the numerical results. 

Another model was calculated assuming an initial composition of 50% amor- 
phous ice and 50% dust by mass and allowing for crystallization of the amor- 
phous ice. The water production rate over the entire nucleus surface is shown 
in Fig. 8 together with the respective temperature map at two heliocentric 
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Fig. 6. Shell temperatures for different depths for comet 67P/C-G with 40° tilt 
angle. The theoretical daily skin-depth for 50% crystalline water and 50% dust 
with Hertz factor of H z = 0.1, was calculated to be 11 cm. The maps indicates that 
the azimuthal temperature homogeneity increases with depth, and that deeper than 
16 cm, hourly temperatures are almost indistinguishable. 

distances post-perihelion. The strong temperature dependence of the subli- 
mation rate is illustrated by the high concentration of the active spot. We 
note that the peak is slightly shifted beyond the subsolar point (toward af- 
ternoon). The maximal integrated total production rate obtained (near per- 
i helion), 8.3 ■ 10 28 s~* ] , is h igher than the observed value of ~ 1 - 10 28 s _1 



([Schleicher and Millisl . 120031 ); this is not surprising since the model assumes 
water sublimation from the entire surface, while active areas are most prob- 
ably confined to a fraction of the nucle us surface, as indic a ted by close-up 



obser vations of comet nucleus surfaces (jKeller et al. 
2006m . 



1986; Sunshine et al. 



4-2 The effect of the spin-axis tilt 



Variations of an input parameter may lead to different behavior patterns for 
models that otherwise have the same properties. Here we show, for the first 
time, the effect of the spin-axis tilt on the temperature distribution over the 
surface of a short-period comet. For illustration, we have chosen the orbital 
parameters of comet 67P/ C-G. Our aim is to determine the minimal and 
maximal temperatures that can be obtained locally, under all possible inclina- 
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Fig. 7. Orbital skin-depth for comet 67P/C-G with 40° tilt angle. For reference, the 
analytical value for 50% crystalline water and 50% dust with H z = 0.1 is 8 m. Top 
maps show a spherical shell, 8 m deep, at perihelion (left) and at aphelion (right). 
Differences between them indicate that this depth is within the orbital skin-depth. 
The bottom maps show temperatures of a deeper shell (12.5 m), where differences 
are minute, indicating that the shell is deeper than the orbital skin-depth. 

tions. This will shed light on the viability of very volatile ices on the surface 
of comet nuclei. 

Assuming a spherical comet nucleus with the orbital parameters of comet 
67P/C-G, we carried out evolutionary calculations over the entire range of 
tilt angles, at intervals of 10°. For all models, the summer solstice point was 
taken to be at perihelion, (i.e. the projection of the spin vector on the ecliptic 
plane, points to the sun at perihelion). Thermal evolution was calculated over 
a period of 5 orbital revolution. The nucleus composition was taken to be 50% 
dust and 50% crystalline H2O ice by mass. Since we were interested in the 
temperature distribution at the surface, and since for a short-period comet, 
amorphous ice is present only at relatively large depths, we have assumed 
crystalline ice throughout. We adopted a 0.1 Hertz factor. 

The results are summarized in Fig. 9. For maximum temperatures (left), south- 
pole temperatures were always lower than the north-pole ones; however, for 
small tilt angles (below 40°) the difference was found to be more significant. 
When the tilted nucleus is at perihelion, the south-pole is directed away from 
the sun, not receiving any solar flux. The only heat flowing into the south 
pole is from non-lateral conduction, which is negligible relative to the solar 
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Fig. 8. Surface temperature and water production rate maps for comet 67P/C-G 
with 45° tilt angle, where crystallization is taken into account. Top maps show 
temperature (left) and local flux (m~ 2 s" 1 ) of water molecules (right) at 1.5 AU 
post-perihelion, and bottom maps, at 3 AU post-perihelion. 



flux. The south-pole remains hidden from solar flux until the nucleus nears the 
equinox, where the sun is directly above the nucleus' equator, and the south 
pole is exposed to the sun at a very low angle. Passing the equinox point, the 
south pole is exposed to direct flux. As the comet retreats from equinox, the 
flux received by the ever exposed south-pole depends on the tilt angle, making 
the small tilt angle pole receive only a fraction of the flux received by a high 
tilt pole. The lowest value of the maximal surface temperature is less than 
130°K. At such temperatures crystallization of amorphous water ice proceeds 
very slowly: the characteristic timescale is of the order of days, competing with 
the dynamical (orbital) timescale. It is thus possible that amorphous ice be 
retained at (or very nearly below) the surface of a very confined area, provided 
that the tilt angle remains unchanged. On the other hand, this temperature 
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Fig. 9. Maximum and minimum local surface temperatures for 10 different spin-axis 
tilts. Extremum temperatures vary only with co-latitude, and are greatly affected 
by tilt angle. The lowest value of the maximum temperature is obtained at the 
southern pole for small tilt angle. The lowest temperatures (over the 5 th revolution) 
are found at the north pole, where the highest maximum temperature occurs. 

is far too high for ices of volatile substances (certainly CO, CH 4 , N 2 , but 
also CO2, NH 3 , HCN) to survive insolation at the surface. These molecules - 
observed to be ejected by comet nuclei — must, therefore originate in deeper 
layers or be released by crystallization of gas-laden amorphous water ice. It 
is, therefore, inconceivable that pristine material be found on the surface of 
comets of the Jupiter family type. 

Minimum temperatures (right panel of Fig. 9) are lowest for the north pole, 
which is allowed to cool for most of the comet's orbit. Temperatures drop below 
60°K, which may render the surface and the outer layers almost completely 
inert. 

To further study the influence of small tilt angles on sub-surface crystalliza- 
tion, a high-resolution model of comet 67P/C-G was computed, adopting a 
15° tilt angle. This high-resolution run of the model took 3 months of calcula- 
tions on a 3.0Ghz 64bit Pentium-4 based PC. Therefore, it was not conducted 
exhaustively for each tilt angle. The model started with 50% amorphous ice 
that was allowed to crystallize. Fig. 10 shows the crystallization cross-section 
pattern (top) and the temperatures for the same area (middle) at 1.5 AU post- 
perihelion. It is shown that the south-pole crystallization front did not exceed 
depth of 3m (during the entire run of the model the front did not advance 
beyond 4m), in contrast to a model with 40° tilt angle (not shown here) that 
reached a depth of 10 m within the same period. The north-pole crystalliza- 
tion front did not exceed the depth of 2 m in both models, although maximum 
temperatures were the highest. This happened due to the rapid exposure of 
the pole to the solar flux. When the pole is hidden from the sun (most of 
the time), the cooling boundary conditions with only small amount of heat 
that has been absorbed during solar exposure, prevent the deeper temperature 
from increasing, and thus crystallization propagates only to a shallow depth. 
This may indicate that pristine material may be found at the north-pole at 
small depths but not at the surface. 
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Fig. 10. Surface temperature and advance of the crystallization front for comet 
67P/C-G with 15° tilt angle at 1.5 AU post-perihelion. The model was started 
with 50% dust and 50% amorphous ice with H z = 0.01 and the amorphous ice 
was allowed to crystallize. Top panels show the amorphous ice mass fraction and 
temperature as function of co-latitude angle and depth, with arbitrary azimuth. 
The crossing white line represents the equatorial plane and the white dot at the 
ordinate represents the sub-solar co-latitude angle. Surface temperature (bottom) is 
shown with same equatorial and subsolar markings. 



We note that our present model does not account for possible recession of the 
surface due to sublimation of the ice. This effect would distort the sphericity 
of the model (numerical grid). On the other hand, as stated above, most 
of the surface of comet nuclei is covered by a dust mantle, presumably of 
high porosity, and water vapor is produced and ejected from an ice-rich layer 
lying beneath the dust mantle. Therefore, the process of erosion may be more 
complicated than just simple ablation and requires separate investigation. The 
point we wish to stress in the present calculation is that amorphous ice may 
lie very close to the nucleus surface under appropriate conditions, although it 
would not be detectable on the surface itself. 
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5 Summary and Conclusions 



We have developed a fully 3-D thermal evolution code for comet nuclei that 
may, in fact, be applied to any small body of the solar system, small enough 
for self-gravity to be negligible. In the present paper, which is the first of a 
series, we have tested the code with respect to analytical solutions, as well as 
1-D and quasi-3-D comet nucleus models, and found excellent agreement. 

We computed models of comet 67P/C-G, adopting a 1:1 mass ratio of ice 
to dust. In this case, surface temperatures are mainly determined by the ice, 
through the strongly temperature-dependent sublimation term of the bound- 
ary condition. Assuming a spin axis tilt of 40°, we obtained a maximum surface 
temperature of 205°K and a minimum temperature of 54°K. 

We have investigated the effect of the spin axis tilt on the surface temperature 
distribution and found that conditions for preservation of pristine amorphous 
ice and moderately volatile species at shallow depths is possible for low tilt 
angles, when some surface areas are shielded from insolation. This is despite 
lateral heat conduction driven by the strong temperature variations at the 
nucleus surface. Although not shown here, more distant comets with low tilt 
angles could retain unprocessed surface amorphous ice. 

A fully 3-D comet nucleus model opens a wide range of possible subjects of 
study, involving inhomogeneities, both innate and evolutionary, which appear 
to be so common to comets. We shall address some of them in upcoming 
papers. 
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